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ABSTRACT 

We quantify quasar color- variability using an unprecedented variability database - ugriz photometry 
of 9093 quasars from SDSS Stripe 82, observed over 8 years at '^60 epochs each. We confirm previous 
reports that quasars become bluer when brightening. We find a redshift dependence of this blueing 
in a given set of bands (e.g. g and r), but show that it is the result of the fiux contribution from less- 
variable or delayed emission lines in the different SDSS bands at different redshifts. After correcting 
for this effect, quasar color- variability is remarkably uniform, and independent not only of redshift, 
but also of quasar luminosity and black hole mass. The color variations of individual quasars, as they 
vary in brightness on year timescales, are much more pronounced than the ranges in color seen in 
samples of quasars across many orders of magnitude in luminosity. This indicates distinct physical 
mechanisms behind quasar variability and the observed range of quasar luminosities at a given black 
hole mass - quasar variations cannot be explained by changes in the mean accretion rate. We do find 
some dependence of the color variability on the characteristics of the fiux variations themselves, with 
fast, low- amplitude, brightness variations producing more color variability. The observed behavior 
could arise if quasar variability results from fiares or ephemeral hot spots in an accretion disc. 
Subject headings: quasars: general - quasars: emission lines - galaxies: nuclei - galaxies: active - 
accretion, accretion discs 



1. INTRODUCTION 

Quasars, the brief phases of high accretion onto the 
massive black holes in the centers of large galaxies, have 
proven to be one of the most versatile classes of astro- 
physical objects in the exploration of the distant Uni- 
verse. Through large efforts and dedicated searches (e.g., 
Schmidt & Green 1983; Croom et al. 2001; Eyer 2002; 
Richards et al. 2002, 2004; Atlee & Gould 2007; Richards 
et al. 2009; D'Abrusco et al. 2009; Bovy et al. 2011a) 
large samples of quasars are known today and have been 
explored in much detail to aid the understanding in fields 
as different as mass clustering on both large and small 
scales (Croom et al. 2005, 2009; Shen et al. 2007, 2009; 
Ross et al. 2009), the understanding of the molecular 
gas content in distant galaxies (Yun et al. 1997; Riechers 
2007a, b), and estimates of cosmological parameters and 
the dark energy equation of state (e.g., Scranton et al. 
2005; Giannantonio et al. 2008; Xia et al. 2009). 

Much effort has been put into understanding the na- 
ture of the quasars themselves and active galactic nuclei 
(AGN). Among the phenomena to be explained is the 
ubiquitous time-variability of the quasar emission. Sev- 
eral physical processes have been invoked to explain the 
variability of the observed optical emission. Foremost are 
accretion disc instabilities (e.g., Rees 1984; Kawaguchi et 
al. 1998; Pereyra et al. 2006), but also large-scale changes 
in the amount of in- falling material may be important 
(e.g., Hopkins et al. 2006, and references therein). Vari- 
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ous stochastic processes have also been suggested as pos- 
sible causes with less success though. 

There are different approaches on how to sort out 
which variability mechanisms are prevalent under what 
circumstance. The foremost diagnostic is the temporal 
behavior of fiux variations: Each of the mechanisms in- 
duces variability on different timescales, from weeks for 
changes on thermal timescales in the accretion disc, over 
months for superpositions of stochastic processes to sev- 
eral years for viscous changes in the large-scale structures 
of the accretion disc and for lens crossing times. These 
'physical' timescales of AGNs (e.g., Webb & Malkan 
2000; Collier & Peterson 2001, and references herein) can 
be compared to the observed AGN variability timescales. 
The observed variability time-scales span the range from 
a few hours (Stalin et al. 2004; Gupta et al. 2005), pos- 
sibly the result of processes in a jet (Kelly et al. 2009), 
to months and years where quasars are known to typ- 
ically vary > 10% (e.g., Giveon et al. 1999; Colher & 
Peterson 2001; Vanden Berk et al. 2004; Rengstorf et al. 
2004; Sesar et al. 2007; Bramich et al. 2008; Wilhite et 
al. 2008; Bauer et al. 2009; Kelly et al. 2009; Kozlowski 
et al. 2010). However, it is a challenging task to disen- 
tangle the processes based on variability-timescales alone 
to get a clear view of the underlying processes. More- 
over, as AGN generally have power-law-shaped structure 
functions for the temporal variations, i.e., with no par- 
ticular timescales, and because measurements of 'charac- 
teristic variability timescales' are likely dominated, or at 
least infiuenced by window functions, interpretations of 
variability time-scales are always somewhat problematic. 
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Nevertheless, there seems to be a general consensus in the 
community, that the most probable scenario for the ma- 
jority of the observed variability behavior is changes in 
the accretion discs. 

The weekly to yearly variability of quasars (AGN) has 
been exploited for several independent purposes, e.g., for 
reverberation mapping (Petersc ' Kaspi et al. 2005, 
2007) to estimate Eddington ratios and black hole masses 
(Peterson et al. 1998, 2004; Kaspi et al. 2000), and for 
quasar identification (Scholz et al. 1997; Eyer 2002; Geha 
et al. 2003; Sumi et al. 2005). Recently, new algorithms 
for identifying quasars via their variability have been es- 
tablished (e.g., Kelly et al. 2009; Schmidt et al. 2010; 
Palanque-Delabrouille et al. 2010; MacLeod et al. 2010, 
2011; Butler & Bloom 2011; Kim et al. 2011). These ap- 
proaches provide an alternative to color selection alone 
for discovery of quasars (e.g., Richards et al. 2002) in 
large-area multi-epoch surveys like the Panoramic Sur- 
vey Telescope & Rapid Response System (Pan- STARRS; 
Kaiser et al. 2002) and the future Large Synoptic Sur- 
vey Telescope (LSST; Ivezic et al. 2008; Abell et al. 
2009) where spectroscopic confirmation of the millions 
of quasar candidates to be found is unfeasible. 

Multi- wavelength AGN variability data to date provide 
clear albeit somewhat qualitative evidence that quasars 
tend to get bluer when they get brighter (e.g., Giveon et 
al. 1999; Trevese et al. 2001; Trevese & Vagnetti 2002; 
Geha et al. 2003; Vanden Berk et al. 2004; Wilhite et 
al. 2005; Sakata et al. 2011). One practical caveat to 
such conclusions is that almost all photometric 'brighter 
makes bluer' claims are based on fitting in fiux versus 
color, without accounting for the color-magnitude error 
correlation in the modeling; such fitting may lead to 
spurious, or at least biased, results as we show in the 
present paper. To avoid these error correlations such 
analyses should be performed on fiux-fiux space. Here 
we develop a better unbiased fitting procedure to quan- 
tify whether brighter-makes-bluer on short time-scales 
(see also Sakata et al. 2011). 

In the present paper, we carry out a comprehensive 
study of 'color variability' in quasars, i.e., we study how 
flux variability is linked to changes of the (observed) op- 
tical colors. We provide both a detailed empirical de- 
scription of the observed variability, and work on linking 
it to the physics of the central engine by quantifying the 
correlation of color variability with the redshift of the 
quasars, their Mbh and L/I/Edd, and to the temporal 
behavior of the fiux variations. While the temporal be- 
havior has been studied extensively, there are no compre- 
hensive studies of color variability in large quasar samples 
with many epochs of data, which offers great potential 
as a diagnostic of accretion disc physics. 

The Sloan Digital Sky Survey (Stoughton et aL 2002; 
Gunn et al. 2006) Stripe 82 provides an unprecedented 
data base for the present color variability study. It has a 
well defined quasar classification and several epochs (on 
average 60 epochs over 8 years) of precise photometry in 
five bands. 

The paper is structured as follows: after introducing 
the superb data set of Stripe 82 in Section 2, we de- 
scribe the development of an unbiased fitting procedure 
in Section 3. In Section 4 we present the color variabil- 
ity results from applying this procedure to the Stripe 82 
data. In Section 5 we discuss our findings and describe 



how the color variability depends on the light curve vari- 
ability properties. Furthermore, we check for L/LEdd as 
well as Mbh dependences, and compare our results to 
recent accretion disc models before we sum-up and con- 
clude in Section 6. The central findings of the paper are 
reflected in Figures 4, 7, and 9. 

2. SDSS STRIPE 82 DATA 

The Sloan Digital Sky Survey's (SDSS's) Stripe 82 
(S82) is an equatorial stripe 2.5 degrees wide and about 
120 degrees long which has been observed many times 
in the 5 SDSS bands over more than 8 years. The S82 
data base hence provides an unprecedented collection of 
data for variability studies in general (e.g., Ivezic et al. 
2007; Bramich et a^ ^^^^; "-sar et al. 2007, 2010; Bhatti 
et al. 2010) and quasars in particular (Kelly et al. 2009; 
Schmidt et al. 2010; MacLeod et al. 2010, 2011; Butler 
& Bloom 2011). The analysis presented here is done on 
the ^9,000 spectroscopically confirmed quasars in S82. 
They have been selected from the SDSS data archive^ 
as described in Schmidt et al. (2010). Thus, we have a 
sample of quasars with on average 60 observations spread 
over a period of roughly 8 years. To obtain further in- 
formation on each individual quasar, such as for instance 
estimates of the bolometric luminosity (I/boi) and black 
hole mass of the central engine (Mbh), we cross-matched 
our list of objects with the catalog of quasar properties 
presented in Shen et al. (2010). We found a total of 9093 
matches which constitute the catalog we will use in the 
remainder of this paper, unless noted otherwise. This 
corresponds to basically the complete catalog of spec- 
troscopically confirmed quasars in S82, hence, matching 
to the Shen et al. (2010) catalog did not cut down the 
sample much. 

3. FITTING COLOR VARIABILITY IN MAGNITUDE SPACE 

In principle optical color variability of quasars, i.e., 
the tendency of changing color, generally becoming bluer 
when they brighten has been well established (Giveon et 
al. 1999; Vanden Berk et al. 2004; Wilhite et al. 2005). 
However, this color variability has been established and 
quantified by fitting data in color-magnitude space, for 
instance in g versus {g — r) space, which suff'ers from 
co-variances between the color and the magnitude un- 
certainties that have not been accounted for in the past 
analyses. As we describe in Appendix A we have found 
that this may lead to severe overestimates of the color 
variability, especially as the photometric errors are not 
negligible compared to the intrinsic variability ampli- 
tudes. To remedy these biases we fit the color variabil- 
ity in magnitude-magnitude space, and then 'translate' 
them into color-magnitude relations. 

The S82 data presented in Section 2 are unprecedented 
in their combination of time coverage, number of epochs, 
filter bands and sample size, which allows us to take the 
color variability analysis to the next level. Since the 
quasar variability typically shows modest amplitude (a 
few tenths of a magnitude or less) it has been charac- 
terized in previous work by a linear relation in fiux-fiux 
space (e.g Choloniewski 1981; Winkler 1997; Suganuma 
et al. 2006; Sakata et al. 2010, 2011). We make the 
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ansatz that the photometric measurements of each in- 
dividual quasar can be represented by a linear relation 
in ^r-space (and m-space). The gr and ui spaces were 
chosen for high S/N and broad spectral range respec- 
tively (the uz-spdice being too uncertain due to 2;-band 
measurement errors). By calculating the Pearson cor- 
relation coefficient (PCC) for each individual quasar we 
verified that this approach is indeed sensible. Averaged 
over the full sample in ^r-space the PCC = 0.8. A PCC 
of 1 indicates a linear correlation with basically no scat- 
ter. Including an extra parameter in the linear relation to 
resemble the intrinsic scatter when determining Sgr also 
shows that the scatter of the assumed linear relation is 
insignificant. Hence, the data does indeed resemble a 
linear relation quite closely. 

When fitting such data a number of factors need to be 
taken into account: there are comparable errors along 
both axes, e.g., for g and r magnitudes, the errors vary 
widely among different data points, and there are 'out- 
liers' (see Schmidt et al. 2010, for physical, observational 
or data processing reasons). In order to take these fac- 
tors properly into account we have used the linear fitting 
approach, including outlier pruning, laid out in Hogg et 
al. (2010) page 29ff. In practice we identify the set of 
relations that make the data likely outcomes 

r - (r) = s'g,{g - (g)) + b (1) 

by a Metropolis-Hastings Markov chain Monte Carlo 
(MCMC) approach. Thus, we are determining the color 
variability (slope) s^^ and the offset on the r-axis b (mov- 
ing each object to its mean g and r value to improve the 
determination of b) by sampling the parameter space via 
an MCMC chain. From simple algebraic manipulations 
of Equation 1 we have that 

r = s'g,g + b' (2) 
g-r = -{s'g,-l){g-{g))+B (3) 

where b' = (r) — Sg^{g) + b and B = —b -\- {{g) — (r)) are 
constants. These equations give the 'transformation' of 
the fit between magnitude-magnitude space and color- 
magnitude space. If {Sg^ — 1) < 0, i.e., if Sg^ < 1 the 
quasar gets bluer as it brightens. In the remainder of 
this paper we will use 

Sgr = {S'gr - 1) (4) 

as our definition of the color variability. In more general 

terms this corresponds to sx^x^ = — 1 where the 

As refer to the photometric bands. This expression has 
the intuitive interpretation that Sgr = means no color 
variability, i.e., brightness variation at constant color, 
whereas Sgr < accounts for the bluer equals brighter 
trend (the most common trend in the data) and Sgr > 
implies objects that become redder when they brighten. 

In Figure 1 an example of such an MCMC fit to con- 
strain the color variability is shown for the quasar SDSS 
J2037-0112 in gr space. The filled circles represent the 
individual photometric epochs from S82, with the ellipses 
indicating the photometric errors in g and r. The 'best' 
fit from Equation 1 is shown as the blue solid line. The 
blue dashed lines show the 68% confidence interval given 
by the 16th-84th inter-percentile range of the MCMC 
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Fig. 1. — Example of MCMC fitting for the gr color variability, 
given a set of independently measured g and r light curve points 
and their uncertainties. These photometric data for the quasar 
SDSS J2037-0112 are shown as black circles with error ellipses in- 
dicating the photometric errors on the measurements. The solid 
blue line indicates the 'best' MCMC fit (i.e., the parameters for 
which the data are most likely) of the linear relation in Equation 1 
to the data. The dashed blue lines show the 68% confidence in- 
terval of the MCMC. The empty data points with dashed error 
ellipses have a posterior probability of being outliers to the rela- 
tion of > 50%, as described in Section 3; see also Section 3 in Hogg 
et al. (2010). This fitting accounts for the independence of the g 
and r measurements, the widely varying error bars, and 'outliers'. 
In the upper left corner we note the mass of the central black hole 
and the bolometric luminosity taken from the value-added quasar 
catalog presented in Shen et al. (2010) and the spectroscopic SDSS 
redshift of the object. In this manner the gr and ui color vari- 
ability of each of the 9093 Stripe 82 quasars in the sample were 
determined. 

'cloud' of possible fits. The described fitting procedure 
also allows estimating the probability that a given data 
point is an outlier to the obtained relation, i.e., an es- 
timate of the posterior probability that each individual 
observation 'belongs' to the obtained relation (see Sec- 
tion 3 in Hogg et al. 2010). Such outliers can be due to 
for instance weather, bad calibration, image defects etc. 
The data points represented by the open circles with the 
dashed error ellipses in Figure 1 have a posterior proba- 
bility of being outliers to the shown MCMC fit which is 
larger than 50%. On average 8% and 19% of the observed 
epochs were counted as outliers to the obtained relations 
when fitting in gr and ixz-space respectively. Similar fits 
in gr (and ui) space were performed for all 9093 quasars 
in the S82 sample, each resulting in estimates for the gr 
and ui color variability for each object. Note that the 
temporal ordering of the flux points plays no role in this 
analysis. 

4. RESULTS 

In the following subsections we will present the results 
from the investigation of the color variability, s^^, s^i^ 
for the 9093 spectroscopic S82 quasars. 
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4.1. Color Variability in gr 

Figure 2 shows the immediate result of the MCMC 
fitting procedure; the directly observed gr color variabil- 
ity for the 9093 spectroscopically confirmed quasars from 
S82 as a function of their spectroscopic redshifts (left 
panel). It is clear that the vast majority of the quasars 
show color variability Sgr < (i.e., they get bluer when 
they brighten) represented by the shaded region. The 
right panel shows the analogous plot for the ui color vari- 
ability, which we consider further in Section 4.3. 

Figure 2 shows that there is a very significant red- 
shift dependence on the mean observed color variability, 
{sgr/ui){^)' Spectra (Wilhite et al. 2005) and other infor- 
mation suggest that the {sgr/ui){z) behavior seen in Fig- 
ure 2 arises from a general trend of bluer continuum color 
in higher flux states modified by the redshift-dependent 
influence of emission lines in a given observed bandpass. 
As we will show in detail below, such a description is 
consistent with trends in the S82 data. 

The colors of quasars in general are known to have a 
pronounced redshift dependence resembling that seen in 
Figure 2 because emission lines and the continuum af- 
fect static or single epoch colors (e.g., Richards et al. 

2001; Wilhite et al. 2005; Wu & Jia 20. o; .ger et 

al. 2011). For instance, the strong drop in {sgr){z) at 
z ~ 0.95 in Figure 2 corresponds exactly to the redshift 
where the Mgll line moves from the g to the r band. This 
is illustrated in Figure 3, where the left panel shows the 
Vanden Berk et a ( 1) composite quasar spectrum 
with the 5 SDSS bands as they would be positioned if 
the quasar was at z = 0.95. The right panel shows the 
ratio between the emission line flux (Fune; the compos- 
ite spectrum minus the estimated continuum flux) and 
the estimated continuum flux (Fcont; modeled as a sim- 
ple power-law) as a function of redshift. The shift of 
the Mgll line from the g (green) to the r (yellow) band 
is marked. Likewise the dips and bumps in {sgr){z) at 
z ^ 1.85,2.8 and 3.5 in the left panel of Figure 2 are 
attributable to the GUI], CIV and Lya hues shifting be- 
tween the g and r bands respectively. In the right panel 
of Figure 3 only 2: < 2.4 is shown since this is the region 
where a simple power-law approximation of the contin- 
uum is valid. For higher redshift the g and r bands move 
blue- ward of the Lya line. 

In order to isolate the continuum color variability for 
comparison with other quantities such as L/LEdd and 
^BH (Section 5.1), we need to eliminate the source red- 
shift dependence induced by the emission lines. This is 
done by 'emission line correcting' the individual values 
of Sgr and Sui by the quantity 

{s)-{su){z) . (5) 

Here the first term is the mean color variability of -0.17 
(-0.46) in gr {ui) which is our stand-in for the emission 
line free color variability (dashed line(s) in Figure 2). 
The second term is the mean color variability for the 
sample capturing the mean redshift dependence of the 
sample as depicted by the red line in Figure 2 for each 
of the individual quasars, k. 

4.2. Reproducing the Color Variability Redshift 
Dependence with Simple Variability Model 



We now quantify to which extent a simple spectral vari- 
ability model can reproduce the observed redshift trends 
in {Sgr/ui){z)- We do this by integrating a time- varying 
sequence of mock spectra created from the Vanden Berk 
et al. (2001) composite quasar spectrum over the SDSS g 
an r filters as illustrated in the left panel of Figure 3. Af- 
ter decomposing the Vanden Berk et al. (2001) spectrum 
in a continuum and line component, by subtracting the 
estimated power-law continuum from Vanden Berk et al. 
(2001), we varied both the continuum and the lines to 
create a mock time sequence of spectra for which we ob- 
tained g and r light curves and then s^^. By changing 
the slope of the continuum (with a pivot-point in the IR 
to ensure Sgr < 0) and scaling the line response by a 
given amount, a sequence of spectra could be created to 
simulate a variable quasar. The line response was charac- 
terized by the ratio between the total integrated change 
in continuum flux and the total change in line flux over 
the modeled wavelength range 



and could be set free (both lines and continuum can vary 
freely) or be fixed. Several setups for creating the se- 
quence of variable spectra were inspected. Among those 
setups were fixed line contribution with changing con- 
tinuum slope and both continuum and lines changing 
in various ways. For given a the emission lines are as- 
sumed to respond instantly to the continuum variation; 
i.e., in this simplistic approach we ignore any reverber- 
ation time-delay between the continuum and the lines. 
An exploration of this effect to carry out reverberation 
mapping (e.g., Peterson et al. 2004; Kaspi et al. 2005, 
2007; Chelouche & Eliran 2011) using the broad-band 
light curves seems promising in light of Figure 2, but 
is beyond the scope of this paper. Details on the simple 
spectral variability models are given in Appendix B. 

The predictions of the spectral variability models are 
shown in Figure 4 together with the estimated values of 
Sgr^ shaded regions, and the mean redshift dependence, 
{sgr){z), from the left panel of Figure 2. This Figure 
shows that {sgr){z) can be best matched if the (implicitly 
instantaneous) line response is very sub-linear: a = 0.1 
(purple line in Figure 4) is a much better fit than the 
model with a = 1 (red line in Figure 4). It is seen 
that for emission lines that vary in lockstep with the 
continuum by a > 25% the redshift features in {sgr){z) 
are 'inverted'. Actually, unresponsive line fluxes (i.e., 
a = 0), lead to the best match in this model context 
(black dashed curve in Figure 4). Overall, Figure 4 tells 
us that the redshift dependence of the gr color variabil- 
ity is nicely reproduced by a simple spectral variability 
model where the continuum of the spectrum is hardened, 
i.e., its power-law slope is changed so brighter makes 
bluer, and the emission line fluxes in the two bands are 
(instantaneously) unresponsive. We know from detailed 
reverberation studies that emission lines do respond (on 
1/2 year timescales Kaspi et al. 2005, 2007) for quasars 
at z ~ 1. The explanation for a ~ may be that the 
lags in the lines are long enough to introduce a phase 
offset whereby the line is sometimes stronger, sometimes 
weaker than predicted from a tight correlation and in the 
net the correlation gets lost. 
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Fig. 2. — The color variability, SA1A25 9^ (left) and m-space (right) for the full sample of 9093 Stripe 82 quasars (dark gray dots), as 
a function of redshift. The black rectangles show the sample mean and its uncertainties in redshift bins of Az = 0.125, containing at least 
5 data points. The width (a) of the color variability distribution is indicated by the thin error bars. The red curve gives the interpolated 
mean redshift trend, {sgr){z), and the black dashed line indicates the sample median color variability of -0.17 (-0.46) in gr and ui. The 
shaded region indicates Sg^^^^ < 0, i.e., where bluer means brighter. The complex redshift-dependence of color- variability is due to the 
influence of various emission lines (see discussion in Section 4.1 and 4.2). 
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Fig. 3. — The role of emission lines in the color variability. On the left, the Vanden Berk et al. (2001) composite SDSS spectrum is shown 
with the 5 SDSS fllter's response curves (ugriz from left to right; arbitrarily scaled for visibility) as they would fall if the quasar was at 
z = 0.95. In the right panel the expected ratio between the line flux -Fune cind the continuum flux -Fcont for the same spectrum is shown 
with the SDSS g (green) and r (yellow) bands shown as solid curves and the u (blue), i (orange) and z (red) bands shown as dashed curves. 
The redshift at which the Mgll and GUI] lines move from the g to the r band {z ~ 0.95 and z ~ 1.85 respectively) has been indicated. 
This illustrates the redshift dependence of the color variability seen in Figure 2 and described in the text. The redshift range in the right 
panel has been set equal to the one used in Figure 2 to ease comparison. 
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Fig. 4. — Comparison of the observed redshift dependence of Sgr 
with our simple spectral variability models. The mean {sgr)(z) is 
shown by the solid black line (cf. Figure 2). The rainbow colored 
curves show the color variability models with fixed ratios a (see 
Equation 6) between the line variability and continuum variability 
from a ratio of 0.1 (purple) to a ratio of a = 1.0 (red). The vari- 
ability model, where only the continuum varies {a = 0) is shown 
as the dashed black line, providing the best match to the observa- 
tions. Reverberation mapping time delays, (~ 0.3...1 year x (l + z) 
in quasars, e.g., Kaspi et al. 2007) are not taken into account. 

Improved modeling, including a broken power-law con- 
tinuum, explicit treatment of line reverberation and lack 
of variability at higher rest wavelengths as shown in Van- 
den Berk et al. (2004) and Wilhite et al. (2005), would 
be fruitful to carry out, but is beyond the scope of the 
present paper. 

4.3. Color Variability in ui 

The SDSS S82 data offer the opportunity to extend 
this analysis beyond the relatively short spectral range 
covered by g and r, 4770 A to 6231 A in the observed 
frame. We do so by exploring the color variability in 
the u versus i magnitude-magnitude space, which cov- 
ers a spectral range from 3543A to 7625A. We chose to 
use the i-band instead of z because of the significantly 
smaller photometric uncertainties in the z-band. The fit- 
ting procedure was exactly analogous to the case of gr 
color variability as described in Section 3. The right 
panel of Figure 2 shows the estimated ui color variabil- 
ity for the S82 quasar sample. Despite the larger scatter 
in Sui at any given redshift we see similar features such as 
a distinct redshift dependence in {sui){z) superimposed 
on quite dramatic overall color variability of -0.46. It 
is clear that the ui color variability is more pronounced 
than the gr color variability; {sui){z) < {sgr){z). This 
holds true for the ensemble properties as well as for indi- 
vidual objects, as illustrated in Figure 5, where we plot 
the emission line corrected (as described in Section 4.1) 
color variability in gr and ui. The fact that Sui < Sgr 
for almost all objects, implies that there is a relatively 
stronger blueing over the ui spectral range than over the 



Fig. 5. — The color- variability in gr vs. ui color variability 
space, after correcting for the line-induced redshift dependence (as 
described in Section 4.1). The gray shaded region shows where 
s < 0, i.e., where bluer means brighter. The dashed line indicates 
Sgr = Sui for reference: the ui color variability is more pronounced 
on average than the gr color variability. Only the 3111 objects of 
the sample with Ssgr < 0.04 are shown. 

smaller gr range. The same conclusion is reached when 
accounting for the difference in the wavelength baselines 
between gr and ui by normalizing Sg^/^i with the corre- 
sponding wavelength ratio. 

The simplistic model for the redshift dependence of the 
gr color variability described in Section 4.2 and shown 
in Figure 4, gives equally good results for {sui){z). 

5. DISCUSSION 

We now proceed to put the color variability into the 
context of other physical parameters that describe the 
quasar phase and the time dependence of variability, ex- 
ploring to which physical processes color variability may 
be linked. 

5.1. Color Variability as a Function of Eddington 
Luminosity and Black Hole Mass 

All 9093 spectroscopically confirmed quasars have 
matches in the quasar catalog presented in Shen et al. 
(2010), of which 99.9% (9088) have an estimate of the 
bolometric luminosity and 84.1% (7615) have an esti- 
mated black hole mass derived from Mgll (see Shen et 
al. 2010, for details). This allows us to normalize the 
luminosity to the Eddington luminosity (i^Edd)- 

If we now plot the emission line corrected (as described 
in Section 4.1) gr color variability against Lboi and Mbh, 
as shown in Figure 6, it is evident that there is no de- 
tectable relation between the color variability Sgr and 
the L/LEdd oi" ^BH- This is also illustrated in Figure 7 
where a 2D histogram of L/LEdd and Mbh, with the bins 
color coded according to the median Sgr^ is shown: across 
the well sampled range in Lboi and Mbh, the median Sgr 
varies by no more than 0.01 as a function of these two 
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variable about its niean value of -0.17. The 2D histogram 
has been smoothed by a 2D gaussian to reflect the un- 
certainty in luminosity and mass, with the full width at 
half maximum of the smoothing kernel (represented by 
the ellipse in the bottom left of Figure 7) FWHM = 
[FWHM(MBH),FWHM(L/LEdd)] = [0.35,0.24] corre- 
sponding to [cr(MBH), cr(^/^Edd)] = [O.lSdex, O.ldex] . 
Plots similar to the ones shown in Figure 6 and 7 for 
the ui color variability show no significant L^dd or Mbh 
dependence, either. In Figure 6 the full sample, i.e., all 
masses and redshifts are shown. Inspecting smaller sub- 
samples in z (and Mbh) space does not change the pic- 
ture. Hence, we find no correlation between the color 
variability in gr (and ui) with L/L^dd oi" Mbh- More 
broadly, this seems to imply that the overall state of the 
quasar (characterized by L/L^dd and Mbh) plays no sig- 
nificant role in determining the color variability. 

5.2. The Color Variability as a Function of the Light 
Curve Variability Characteristics 

In Schmidt et al. (2010) we characterized the r-band 
variability of all the 9093 quasars through a 'structure 
function' with an amplitude parameter A and the light 
curve stochasticity, 7. The structure function variability 
of each individual quasar was modeled by a simple power- 
law 

SF^od(AUs|A7) = ^(^)^ (7) 

with Atobs being the time between the observation 
of two individual photometric epochs and SFmod = 
a/ ((m(ti) — m(t2)Y). The structure function of a peri- 
odically varying object or one varying like white noise 
will have a fiat structure function and hence a small 
power-law exponent 7. Thus a large 7 indicates a secu- 
larly varying object or an object with a random walk like 
variability. The latter has been shown to describe quasar 
variability well in Kelly et al. (2009) and MacLeod et 
al. (2011). The amplitude A corresponds to the aver- 
age variability on a 1 year timescale. In Schmidt et al. 
(2010) all calculations were done in the observed frame 
to mimic quasar identification with no prior information 
such as redshift. However, the spectroscopic redshift of 
each of the 9093 S82 quasars is known and the ampli- 
tude can be corrected for time-dilation. The rest-frame 
variability amplitude A' is defined to be ^4(1 + 2:)^ such 
that 

SF^od(Airest|^',7) = A' y (8) 

where Arrest is now the difference between observations 
in the quasar rest-frame. All quoted A's are estimated 
from the robust r-band measurements as described in 
Schmidt et al. (2010). The variability amplitudes are 
independent of z in agreement with Givcori et al. (1999) 
and the majority of the previous studies listed in their 
Table 1. 

In the following, however, the structure-function pa- 
rameters A and 7 have been obtained somewhat differ- 
ently from Schmidt et al. (2010). Rather than fitting 
the structure function directly to the magnitude differ- 
ences we fit a Gaussian Process model ( 11 and 
Williams 2006) defined by the structure function to the 
magnitudes directly. This properly includes all of the 



correlations between data points. This Gaussian Process 
model consists of an n-dimensional Gaussian distribution 
(for n epochs) with a constant mean m and n by n vari- 
ance matrix V . The elements of this variance matrix are 
given by 

y.. ^ V{\U-tj\) = V{^U^) = \ [SF2.(oo) - SF2.(At,,)] 

(9) 

for data points at epochs ti and tj. Here the structure 
function SF^j is given by 

SF,,- = sj(im{U) - m{ti)f) . (10) 

The photometric-uncertainty variances are added to the 
diagonal elements of V. For the power-law structure 
function we cut off the power-law at 10 years such that 
SF(oo) is finite. As all data span less than 10 years this 
cut-off does not influence the fit. This type of fit is simi- 
lar to the Ornstein-Uhlenbeck process describing quasar 
variability as a damped random walk (e.g., Kozlowski et 
al. 2010; Butler & Bloom 2011; MacLeod et al. 2011). 
For more details, see Bovy et al. (2011b). 

We can now look at the emission line corrected Sgr 
and Sui as a function of A' and 7 for all quasars. Fig- 
ure 8 shows that Sgr/ui seemingly vary both with A' and 
with 7. However, the limit of little variability (small 
A') requires particular care, both because outliers play 
a bigger role and because A' and 7 starts to be de- 
generate (Schmidt et al. 2010). We estimated the gr 
color variability of 500 color selected non-varying FG- 
stars (see Schmidt et al. (2010) for further details) and 
of 483 S82 RR Lyrae stars from Sesar et al. (2010). These 
are over-plotted in the top panel of Figure 8 as the blue 
and red points respectively. As expected the RR Lyrae 
have a well defined color variability, whereas the inferred 
color variability of the non- varying FG-stars span a much 
wider range of Sgr. Interestingly, the majority of the 
non-varying FG-stars have color variability estimates of 
Sgr < like the quasars and the RR Lyrae stars. This 
seems to be caused by the outliers in g being relatively 
larger than the outliers in r, hence affecting the initial 
guess of the MCMC in a bluer-brighter direction. In 
the case of the RR Lyrae the well defined mean color 
variability in gr is expected as RR Lyrae change their 
effective temperature and luminosity during their pulsa- 
tion. By creating a sequence of black body spectra with 
temperatures from 6200K to 7200K, estimating the flux 
received in the g and r bands for each spectrum, and us- 
ing that as a simple model for a variable RR Lyrae star, a 
color variability of Sgr ~ —0.23 is obtained, in very good 
agreement with the observations (Figure 8, top panel, 
red dots). Thus, in general the Sgr for the FG and RR 
Lyrae stars look as expected. 

A more direct way to estimate the fidelity of color vari- 
ability estimates at low A' values is to recover Sgr esti- 
mates for objects of known (simulated) color variability. 
We induce such simulated variability into the 500 FG- 
stars by generating data of a certain 1 year amplitude 
(A^) from the original FG-star g and r photometry. We 
do that by generating new u^g^r and i magnitudes for 
the individual epochs j via the expression 

AMJD, 

Pi,sim — Pj + ^sim Q^r 5 (Hj 
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Fig. 



The emission line corrected gr color variability does not depend on L/L-^dd (left) or on M-qu/Mq (right). The black rectangles 
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indicate the mean Sgr and its error in each L and Mbh bin. The thin error bars show the width of the Sgr distribution cf. Figure 2. No 
correlation between the color variability and L/L-^dd cind M-qh/Mq is detected. This implies that the color variability is not a function of 
black hole mass or the overall accretion (Eddington) rate. 

panels of Figure 8 this recovered mean average (over 50 
randomly chosen FG stars) color variability is shown for 
a sequence of variability amplitudes, A'^i^^ as large filled 
green circles. This shows that the recovered gr (ui) color 
variability has some systematic errors below < 0.1 
{A^ < 0.25). When ignoring quasars with variability am- 
plitudes smaller than 0.1 (0.25) the correlation between 
Sgr and Sui and A' is still present and significant. The 
trustworthy part of the relations, (5^^/^^) (A, 7), in Fig- 
ure 8 is shown in Figure 9. Again the black rectangles 
represent the uncertainty on the estimated mean color 
variability. The left hand plots correspond directly to 
Figure 8, whereas for {sgr/ui){l) ^^e black rectangles are 
estimated only from objects with A^ > 0.1 {A^ > 0.25). 
Figure 9 clearly illustrates the trend that objects with 
larger variability amplitude have a smaller color variabil- 
ity (meaning less blueing when brightening) than for low 
A\ On the other hand the color variability is indepen- 
dent of 7. 

As mentioned Vanden Berk et al. (2004) and Wilhite 
et al. (2005) showed that there is a lack of variability 
at high rest wavelengths. Furthermore, it is known that 
quasar variability is anti-correlated with luminosity (e.g.. 
Hook et al. 1994; Cristiani et al. 1996; Vanden Berk et al. 
2004). This might lead to the suspicion that the trend 
presented in Figure 8 and 9 is nothing more than a red- 
shift effect. If this was the case the relation should be due 
mainly to \ow-z objects, since the most variable quasars 
are supposedly low luminosity quasars, i.e., necessarily 
only observed at low and should therefore disappear 
at high redshift. Estimating the relation between Sgr and 
A^ in various redshift bins (also split in luminosity) shows 
that the trend is equally strong for all redshifts and all 
luminosities. Hence, the presented relation appears to be 
of a physical origin and not merely a redshift effect. 




8.5 9.0 
log(MBH/M^ 

Fig. 7. — The color- variability Sgr as a function of L/L-^dd cind 
5 as in Figure 6, again showing no significant trends. 
The color coding indicates the emission line corrected median 
gr color variability of the objects in the bin. To refiect the 
uncertainty in luminosity and mass, the distribution has been 
smoothed by a 2D gaussian (represented by the ellipse in the bot- 
tom left corner) with full width at half maximum of FWHM = 
[FWHM(MBH),FWHM(L/LEdd)] = [0.35,0.24] corresponding to 
[fT(MBH),fT(L/LEdd)] = [0.15dex,0.1dex]. 

by construction a data set with 5 = 0. Here p repre- 
sents the photometric measurements in a given band, j 
runs over the individual epochs, and AMJDj refers to 
the observation time of the jth epoch with respect to 
the first observation. In this way many of the aspects of 
the real data (i.e., the outliers and realistic photometric 
errors) are included in the simulated data. In the two left 
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Fig. 8. — Correlation between the light curve (temporal variability) structure function and color variability. The emission line corrected 
color variability in gr (top) and ui (bottom) are plotted against A' (left) and 7 (right) from Schmidt et al. (2010). A' indicates the mean 
level of variability within one year (rest-frame), and 7, the power law exponent in the structure function, indicates how random (low 7) or 
secular (high 7) the light curve variations are. Gray dots (data), black rectangles and error bars are analogous to Figure 6. In the gr plots 
(top) 500 non- varying FG stars and the 483 RR Lyrae from Sesar et al. (2010) are shown on top of the 9093 S82 quasars (gray dots) as 
blue and red points respectively. In the left column the recovered average Sgr (sui) for 50 FG stars is shown as green filled circles, where 
synthetic brightness variations with Spr(in) = and different variability amplitudes A^^^^^ had been created. As described in the text this 
illustrates that only gr (ui) trends for > 0.1 (A' > 0.25) can and should be trusted. It shows that objects with large A' , i.e., with large 
variability amplitudes, have a color variability close(r) to 0, i.e., less blueing when brightening, than do objects with small A' . The trends 
in the two right hand plots are dominated by low A' objects and is therefore not trustworthy. 



5.3. Color Variability and Changes in the Mean 
Accretion Rate 

In this Section we carry out a cursory exploration as 
to the physical origin of the observed color- variability. 

5.3.1. Color Variability of Individual Quasars vs. the Color 
Distribution of Quasar Ensembles 

The colors of quasars at a given redshift are known 
to depend only weakly on their mean accretion luminos- 



ity or accretion rate (Davis et al. 2007), while we find 
that individual quasars become considerably bluer when 
they brighten on year time-scales. This suggests differ- 
ent physical mechanism creating the accretion luminosity 
range in ensembles and the luminosity variations in indi- 
vidual quasars. 

In Figure 10 the emission line corrected color variabil- 
ity of a sub-sample of the S82 quasars is shown in gr 
and iii-space. This sub-sample represents the 'average' 
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Fig. 9. — The estimated correlations between mean color variability and structure function parameters. The four panels show the portion 
of parameter space in Figure 8 that is well populated by the sample. The rectangles in the left hand plots correspond to the ones in 
Figure 8. In the right hand plots the rectangles are only estimated from objects with A' > 0.1 (0.25). A clear trend is seen, that the 
color variability of objects with high variability amplitude (left panel) are closer to than objects with small variability amplitude. This 
indicates that strongly varying quasars get less blue as they brighten on average, than do moderate varying quasars. The color variability 
is independent of the power law index 7. 



quasars, i.e., the combined sample of the 33rd-66th per- 
centile of masses and the 25th-75th percentile of redshifts 
for the quasar sample. The color variability of each indi- 
vidual quasar is depicted as a short solid gray line show- 
ing Sg^ W^i) from Equation 1 for each quasar centered on 

[{9)^ (^)] {[{^)^ i^)]) ^^^^ particular quasar. Only every 
10th object of the sub-sample is actually shown to keep 
the individual gray lines visible. The length of the lines 
resembles the change in the photometric ^-band (ix-band) 
data of the quasar. The Figure compares the average Sg^ 
(sJ^J of all the individual quasars in the sub-sample (red 
solid line) with a fit to the time-averaged color distri- 



bution of the sub-sample (black solid line) where each 
data point corresponds to [(^), {r)]k {[{u)^ (^)]/e) with k 
counting the quasars. Figure 10 reveals that indeed the 
mean color variability for individual quasars is much 
more pronounced than the equivalent quantity for the 

ensemble Sgampie = d\ni^l ' given sub-sample 

-0.18 and Sui ^ —0.48 on average (as opposed 



to s 



gr 



-0.17 and s., 



compared to s 



gr 



-0.01 and s., 



-0.46 for the full sample) 



-0.08 for the cor- 



responding time-averaged sub-sample color distribution. 
The difference is highly significant in both cases, with the 
ui color variability difference formally larger, because of 
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the broader spectral range. The exact same trends are 
found for plots containing the full quasar sample. 

This result shows that (temporal) color variability of 
individual quasars is considerably stronger than the color 
range of ensembles of quasars at similar redshifts and 
with similar black hole masses, that presumably differ in 

5.3.2. Color Variability vs. Accretion Disc Models 

Explaining quasar spectral energy distributions, and in 
particular the optical/UV continua through steady-state 
accretion disc models has an established history (e.g., 
Shakura & Sunyaev 1973; Bonning et al. 2007; Davis et 
al. 2007). However, comparing the observed color vari- 
ability of large samples of quasars with the predicted col- 
ors of model sequences of varying accretion rate has not 
been done yet. Such a comparison could tell us whether 
it is sensible to think of the quasar variability on scales 
of years as changes in the mean accretion rate. The su- 
perb S82 data enables us to perform such a comparison, 
by comparing the observed color and color variability of 
the S82 quasars with sequences of accretion disc models 
presented in Davis et al. (2007). 

Davis et al. (2007) presented three different thin ac- 
cretion disc models that describe the spectral slope of 
quasars as a function of i^boi/^Edd and Mbh- We took 
these three models and worked out predictions for the 
observed g and r band for models of a given Mbh but 
varying accretion rates. The three models presented in 
Davis et al. (2007) and the color we adopt for their graph- 
ical representation are: 

1) A relativistic model of accretion onto a 
Schwarzschild black hole with a spin parame- 
ter of 0. The emission is based on Non-LTE 
atmosphere calculations (green). 

2) A relativistic model of accretion onto a 
Schwarzschild black hole with a spin parame- 
ter of 0. The disc is emitting as a black body 
(red). 

3) A Model of accretion onto a spinning black hole 
(spin parameter of 0.9) with emission based on 
Non-LTE atmosphere calculations (orange). 

For further details on the models we refer to Davis et al. 
(2007). 

We can then compare the models to the data in two 
respects: do they predict the right color (which has been 
done before) and do they predict the right change of 
color with changing accretion rate or luminosity? In Fig- 
ure 11 the object from Figure 1 is shown in g-r-{g — r) 
space (without error ellipses) together with its best fit 
color variability (blue solid line). The three accretion 
disc models are shown as solid lines in bundles of three, 
where each of the three lines corresponds to a different 
black hole mass, as noted in the bottom right corner of 
each panel. In this particular case model 2) matches the 
data well both in color and in the change of color with 
changing luminosity. However, such a good match is not 
representative for the ensemble. We quantify this for the 
whole sample by estimating the 'goodness' of the models 




as: 



j, model ^-|^2^ 

^ , , . ' O'Ui 

3 •'I J 

As = Sk- 5/e,model (13) 

where y = {g — r) and x = g. The index j runs over 
the Nk epochs for each individual quasar k. The model 
prediction is the color at a given luminosity (or accre- 
tion rate), for a fixed black hole mass y{x)j^^ode\' The 
photometric error on the color for the jth measurement 
is denoted as 5yj. Since the error on the Mbh estimates 
based on Mgll (Shen et al. 2010) is 0.4dex (Vester- 
gaard & Peterson 2006; De Rosa et al. 2011), we have 
chosen to show three values for Mbh, leading to three 
model prediction lines, for each model in Figure 11. 

The 'goodness' parameters and As defined in equa- 
tions 12 and 13 therefore describe how well the observed 
values of color and color variability are predicted by the 
Davis et al. (2007) models. D^^ can be seen as the stan- 
dard measure of comparison between model and data 
before squaring, i.e., it estimates the difference between 
the model color and the observed color averaged over all 
epochs for each quasar. The As is simply the difference 
between model color variability and observed color vari- 
ability of each quasar. 

Figure 12 summarizes the model data comparison for 
a subset of quasars with z ~ 1.5: the quantities from 
equations 12 and 13. All models predict a color variabil- 
ity - as a function of changes in the mean accretion rate 
- that is weaker than observed. On average model 2) 
shown in red matches the observed gr color variability 
the best. Furthermore, the D^_^ values indicate that 
the g — r color is on average overestimated by model 1) 
and 3), whereas the distribution of model 2) has a mean 
very close to the dashed perfect agreement line. Creating 
similar plots for other mass and redshift ranges as well 
as for the results in i^i-space show the exact same trends. 
Thus, of the three Davis et al. (2007) accretion disc mod- 
els considered here, model 2) matches the observed color 
and the obtained gr and ui color variability the best. 

6. CONCLUSION 

In the present study we determined and analyzed 
the color-variability of 9093 spectroscopically confirmed 
quasars from SDSS Stripe 82, to understand to which ex- 
tent and why quasars get bluer (redder) if they brighten 
(dim), by fitting linear relations between the SDSS g 
and r bands as well as between the u and i bands in 
magnitude-magnitude space. The connection of various 
quasar properties to the color variability were inspected 
before the results were compared to models of accre- 
tion disks with varying accretion rates from Davis et al. 
(2007). Our main results can be summarized as follows: 

1. We showed that quasar color variability, sx^x^ = 

f^^^ — 1, is best determined by fitting data in 

the statistically independent magnitude-magnitude 
space, rather than in color-magnitude space as 
many studies have done. Unless care is taken to 
account for the data correlations, the latter ap- 
proach may lead to spurious or biased estimates of 
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Fig. 10. — The color-variability in individual quasars (red) vs. the color-luminosity relation in the ensemble of quasars (black), drawing 
on the sub-sample of 33rd-66th percentile of masses and the 25th-75th percentile of redshifts. The mass and redshift ranges are shown in 
each panel. Each individual quasar is indicated by a short solid gray line corresponding to the fitted and emission line corrected gr (ui) 
relation, with a length refiecting the standard deviation of the g (u) light curve. Only every 10th object of the sub-sample is shown to keep 
the individual objects visible. The red solid lines show the average trend within this sub-sample of Sgr ~ —0.18 and Sui ~ —0.48 with the 
dashed red lines indicating the bootstrapping uncertainty. The black lines show an MCMC fit to the time-averaged fiuxes for the sample, 
i.e., the fit to ({g), (r)) ± (stdev(gf), stdev(r)) for each quasar corresponding to Sgr ~ —0.01 and Sui ~ —0.08. The dashed black lines show 
the 68% confidence interval of the MCMC fit. In other words, the red line is the ensemble mean color variability while the black line is a fit 
to the ensemble of time-averaged mean magnitudes. Hence, it is clear that the average color variability of the individual quasars deviate 
significantly from the time averaged sample color variability. 




Fig. 11. — Comparison of the observed color variability with sequences of steady-state accretion disc models (Davis et al. 2007) of different 
accretion rates. This comparison is illustrated using the quasar SDSS J2037-0112, also shown in Figure 1, shown in r vs. g (left) and 
g — r vs. g (right) space. The black symbols are the individual photometric measurements. The blue solid line shows the MCMC fit in 
gr-space to these measurements and the blue dashed lines indicate the 68% confidence interval of that fit. The open symbols denote likely 
outliers (see also Figure 1). The Davis et al. (2007) models described in Section 5.3.2 are shown in green (model 1), red (model 2) and 
orange (model 3) with L/LEdd changing along the lines. Each set of 3 lines denotes models for differing black hole masses near the value 
determined from the Mgll line width (Shen et al. 2010). From top to bottom each set of lines (models) correspond to 9.1, 9.0 and 8.9 
log(Mniodel/^0)- For this object the match with model 2 is good; see Figure 12 for an ensemble comparison. 
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Fig. 12. — Comparison of the observed color variability to steady- 
state accretion disc models with changing L/LEdd Each individual 
point reflects the data model discrepancy in the mean color Dg_r 
(Equation 12) and Asgr (Equation 13) for one individual quasar 
in the parameter space illustrated in Figure 11. The plot shows 
the sub-sample of 1066 quasars with z = 1.4-1.6. The color of the 
points correspond to model 1) (green), 2) (red) and 3) (orange) 
from Davis et al. (2007) converted into g and r magnitudes. The 
dashed lines show perfect agreement between model and data in 
color and color variability, respectively. The projections of the 
distributions are shown as histograms along the axes. While the 
different sets of models can reproduce (by design) the mean quasar 
colors (i.e., Dg_r), the observed color variability is far too strong to 
be interpreted as changes in accretion rate in a steady-state model 
context. 

color variability. The gr and ui color variability for 
the vast majority of the 9093 quasars Sg^/ui < 0, 
confirming that quasars get bluer when they get 
brighter. 

2. The color variability as measured in gr and ui space 
exhibits a distinct redshift dependence, which we 
could clearly attribute to the effect of emission 
lines exiting/entering the photometric SDSS bands. 
From a set of simple models of spectral quasar color 
variability, a model in which the line and contin- 
uum vary in phase but with the line amplitude fixed 
to 10% or less of that of the continuum, is able to 
reproduce the observed redshift trends in the gr 
(and ui) color variability as well as the observed 
amount of color variability. 

3. The fact that we see clearly the impact of the 
emission-line fluxes on the broad band photome- 
try through the redshift-dependence of the color 
variability implies that broad-band reverberation 
mapping should be possible with the data set at 
hand. 

4. Correcting for the emission lines leaves us with a 
sample mean (continuum) color variability of Sgr = 
^ - 1 = -0.17 and, analogously, Sui -0.46. 

5. We found that the emission line corrected color 
variability is independent of L/LEdd and Mbh in 



both gr and ui: there is no correlation between 
the mass and luminosity of quasars and their color 
variability. 

6. The color variability, however, does depend on the 
light curve variability properties (described by a 
power-law structure function as in Schmidt et al. 
2010). We found that quasars with large variability 
amplitudes {A') tend to have less color variability, 
as compared to quasars with small variability am- 
plitudes. 

7. We found that the characteristic color variability 
on timescales of years of the individual quasars 
is larger than the dependence of the typical 
quasar colors on their overall accretion state (i.e., 
L/LEdd)- This implies that changes in the overall 
accretion rate cannot explain the observed color 
variability. Ephemeral hot spots may however be 
a plausible explanation for the observed color vari- 
ability. This picture is confirmed by our compar- 
ison of the observed color variability to sequences 
of steady-state accretion disc models by Davis et 
al. (2007) with varying accretion rates, which also 
exhibit much less color variability as a function of 
accretion rate. 

Our analysis provides a clear indication that on time- 
scales of years quasar variability does not reflect changes 
in the mean accretion rate. Some other mechanism must 
be at work; presumably some disc instability. What 
mechanism match the existing data, certainly warrants 
further modeling. The current study can also be viewed 
as an initial foray into the realm of multi-band, multi- 
epoch panoptic photometry that the Pan- STARRS and 
LSST surveys can bring to full fruition. 
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Fig. 13. — Simulated data illustrating the effect of error covariances between magnitude and color, that has not been accounted for in a 
number of previous analyses. The left panel shows 5 'clouds' of simulated observations in r-g space drawn from 2D gaussian distributions 
at = 18, 19, 20, 21 and 22 with 'photometric' errors of 0.02, 0.025, 0.04, 0.06 and 0.15 respectively. In the right panel these data are 
plotted in the g-{g — r) color-magnitude space. The stripy pattern seen here is due to the correlation between the errors on the g magnitude 
and the g — r color. Our test showed that this effect need to be accounted for when analyzing data of SDSS Stripe 82 quality, as is also 
illustrated in Figure 14. 
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APPENDIX 

A. FITTING IN COLOR-MAGNITUDE AND MAGNITUDE-MAGNITUDE SPACE 

As mentioned in the text the photometric errors of the {g — r) color are correlated with the photometric errors of the 
g and r band. Thus, when estimating the color variability of objects in general, and quasars in particular, as is done 
in the present paper, care has to be taken that the co- variances of the errors are either removed, taken into account or 
avoided. Above we have avoided the co-variances by estimating the color relation in magnitude-magnitude space and 
then 'translated' that into a color variability in color-magnitude space as described in the text. In the following we 
illustrate the problems one can run into, if the color variability is instead estimated directly in color-magnitude space. 

The correlated errors between for instance the g band and the {g — r) color are easily illustrated, by simply drawing a 
set of random 'observations' from a gaussian distribution with a standard deviation corresponding to the approximate 
photometric error at the given magnitude. In the left panel of Figure 13 such a sequence of simulated data is shown. 
The data have been drawn from 2D gaussian distributions in g and r with mean magnitudes of approximately 18, 19, 
20, 21, and 22 and estimated errors of 0.02, 0.025, 0.04, 0.06, and 0.15 respectively. The soHd line shows the g = r 
relation for reference. The color trend (deviation of the data from the g = r line) has been put in to mimic the average 
color trend of quasars at the given magnitudes. Plotting these simulated observations in color- magnitude space, as 
done in the right panel of Figure 13, clearly illustrates the error correlations. The stripy pattern of the color-magnitude 
diagram is not a consequence of a color change in the object, but a consequence of the fact that the errors in {g — r) 
are correlated with the errors in g\ this is the reason that the 'length' (or artificial color change) of each set of points 
grows for fainter magnitudes. 

In Figure 14 we show how this effect looks when dealing with real data. Figure 14 shows the SDSS S82 photometric 
data of the quasars SDSS J0320-0051 (left) and SDSS J2141-0050 (right) in magnitude-magnitude, magnitude- color 
and color-magnitude space. The top panel corresponds to the left panel of Figure 13. The stripy nature of the data 
when turned into colors is clearly visible in the center and bottom panels of Figure 14. The data have been color 
coded according to the observation time. In each panel one solid and two dashed lines are shown. The solid line has 
been fitted to the shown data, whereas the dashed lines are 'translations' of the fits from the other two spaces. It is 
clear that in the case of the superb data set of SDSS S82 the difference between fitting in r-g (top), g-{g — r) (center) 
and {g — r)-g (bottom) space is negligible. However, if one imagines that only data from year 6 and 7 were available 
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for SDSS J2141-0050 (right panel) and the color variability was estimated based on either g-{g — r) or {g — r)-g, it is 
clear that the fit would deviate significantly from the 'real' color relation because of the co-variant errors. 

If this effect is not taken into account or avoided when estimating color variability, and the error co- variance 'color 
change' is interpreted as an actual color change of the quasar, there is a high probability that the results and conclusions 
will be erroneous. As mentioned in the text color variability has been estimated in color-magnitude space in the past 
with only a few exceptions. Hence, there might be cases in the literature where this effect has not been taken properly 
into account and therefore might affect the validity of the results. It is hard to quantify how much this effect will 
affect the results and conclusions made so far in the color variability literature, and we have therefore not made any 
attempts at quantifying it, but will just note that one needs to take this error correlation into account or, as it is done 
here, estimate the color variability in magnitude-magnitude space to avoid it. 

B. SIMPLISTIC SPECTRAL VARIABILITY MODEL 

In the following we describe the simplistic model for the spectral (color) variability of the quasars used in Section 4.2. 
The simple spectral variability model is an attempt to reproduce the observed redshift trends in the mean color 
variability, {sgr/ui){^)- The model is based on the composite SDSS spectrum from Vanden Berk et al. (2001), -FvdB- 
We decomposed the composite spectrum in a line component, i^une? and a continuum component, i^cont? such that 

^VdB = ^line + ^cont • (Bl) 

The underlying continuum of the composite spectrum is well modeled (for 1150A < A^est ^ 4500A) by a simple 
power-law with a power-law index of = —1.528. By fixing the power-law continuum model with a pivot point 
in the IR (ensuring that 5^^/^^ < 0) we simulate the variable continuum by a (time) sequence of power-laws with 
different Px. Adding fractions of Fiine to each power-law simulates the (potentially) variable emission lines. The 
amount of variability in the emission lines can be fixed to the variability of the continuum power-law via a constant a 
in Equation 6, which as a reminder reads 

^-^line /^\ 



We define 



^Fiine= / Fiinei^.tj) dX — / Fiine(A, tj_l) C^A 

P^max P^max 

^^cont = / Fcont{X, tj) dX — I 



(B2) 



ont(A,t,-i)^iA (B3) 

where Amin ^ 1150A and A^ax ^ 4500A. The tj and refers to the 'epochs' of the variability model. This implies 
from Equation 6, that for a fixed a the emission line response for each variability model 'epoch' is given by 

P^max P^max P^max 

/ Fiine(A, tj) dX = / Fiine(A, t,_i) dX ^ a Fcont(A, tj) - F^ont{\ tj_i) dX . (B4) 

Here the emission lines are assumed to respond instantly to the continuum variation. By integrating the obtained 
spectra at the different epochs, tj, over the SDSS bands the model color variability can be estimated. The results 
shown in Figure 4 are for a variability model with fixed a, since such a variability seems to resemble the observed 
redshift dependence of the color variability the closest. 

The pivot point for the results shown in Figure 4 was put at 4Amax ^ 1.8/im. Changing the pivot point slightly 
changes the amplitude of the obtained model predictions (dashed and colored curves in Figure 4) and the mean color 
variability in such a way that pivot points in the far-IR result in smaller amplitudes and larger {sgr/ui){^)- The actual 
curvature of the curves in Figure 4 does not change with the pivot point, i.e., the redshift dependence is independent 
of the continuum power-law model pivot point. 

This spectral variability model predicts the redshift dependence in Sgr and Sui equally well. 

A model similar to the one presented here, was used in Richards et al. (2001) to explain the redshift dependences of 
the SDSS quasar colors. 
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